function hypercolumn(p_filename)
% Version infomation:
% 1 - Weighting the matrix in the connectivity function
% 2 - Wichten the selectivity smax unterschiedlich fuer excitatory und inhibitory
% 3 - 09.05.13 Change structure of p 

    %% close and clear all
    close all; 
    warning off;  
    addpath('input');
         
    %% Define global variables
    global  x y nmass te ti he hi e0 r0 v0 n sz neqn tend ...
        c1s c1n c1l c2s c2n c2l f1 f2 f3 f4 f5 p sample_factor path 
     
    %% Load initial data file  
    % te=p(1); ti=p(2); he=p(3); hi=p(4); 
    % sz = p(5); sample_factor = p(6)
    % c1s=p(7); c1n=p(8); c1l=p(9);
    % c2s=p(10); c2n=p(11); c2l=p(12);
    % f1 = p(13); f2 = p(14); f3 = p(15);
    % f4 = p(16); f5 = p(17); tend = p(18);
    % disp('Loading initial parameter...');
    %load(['param/' filename '.mat']);
    load([p_filename '.mat']);
    fprintf('FILE:[%s] C:[%.2f %.2f %.2f %.2f %.2f %.2f] F:[%.1e %.1e %.1e %.1e %.1e]\n',...
        p_filename,p(7),p(8),p(9),p(10),p(11),p(12),p(13),p(14),p(15),p(16),p(17));     
    
    %% Parameters
    % synaptic and dendrite parameter
    te      = p(1);% 10  % ms
    ti      = p(2);% 20; % ms
    he      = p(3);% 3.25;  % mV
    hi      = p(4);% -22;% mV
    sz      = p(5);% integration stepsize = 0.1 ms
    % Asymmetric weighting factors
    c1s     = p(7);% Weighting exc. = 1
    c1n     = p(8);% Weighting exc. next state = 1.25
    c1l     = p(9);% Weighting exc. last state = 0.75
    c2s     = p(10);% Weighting inh. = 1
    c2n     = p(11);% Weighting inh. next state = 0.75
    c2l     = p(12);% Weighting inh. last state = 1.25
    % Connectivity weigting factors
    f1      = p(13);% Scaling factor of Wpp (distal) (not in 1-Hypercolumn-Model)
    f2      = p(14);% Scaling factor of Wpi (distal) (not in 1-Hypercolumn-Model)
    f3      = p(15);% Scaling factor of Wip (local) (not in 1-Hypercolumn-Model)
    f4      = p(16);% Scaling factor of Wii (local) (not in 1-Hypercolumn-Model)
    f5      = p(17);% Scaling factor of Wpp (local) (not in 1-Hypercolumn-Model)
    tend    = p(18);% total integration time ms = 5000
    sample_factor = p(6); % = 0.1 -> maps 9x6; 0.2 -> maps 18x12; 0.1 -> maps 30x20
    % Parameter of sigmoid function Jansen and Rit
    e0       = 2.5e-3;% 1/ms
    r0       = 0.56;  % 1/mV
    v0       = 6;     % mV
    % others
    x        = 1;% number of hypercolumn x direction
    y        = 1;% number of hypercolumn y direction
    n        = x*y;% number hypercolumns
    nmass    = n*3*8;% number JR masses
    neqn     = 2*nmass;% number ODE
       
    % Result path
    path     = 'output';% where to save the results

    %% Solving the ode system
    [OM] = VC_ode();
    %fname = [path '/c1s=' num2str(c1s) '_c1n=' num2str(c1n) '_c1l=' num2str(c1l)...
    %    '_c2s=' num2str(c2s) '_c2n=' num2str(c2n) '_c2l=' num2str(c2l) ...
    %    '_x=' num2str(x) '_y=' num2str(y) ...
    %    '_f1=' num2str(f1) '_f2=' num2str(f2) '_f3=' num2str(f3) '_f4=' num2str(f4) '_f5=' num2str(f5)];
    %save([fname '.mat'],'OM');
    %plotFigure(fname,OM);
    
    load input/cmap8OMs.mat   % color map
    OM = OM(:,end-4000:end);
    for i = 1:8
       subplot 411; hold on; plot(OM(i,:),'color',cmap(i,:)/255); axis tight; title('PSP Pyramidal');
       subplot 412; hold on; plot(OM(i+8,:),'color',cmap(i,:)/255);axis tight; title('PSP Excitatory');
       subplot 413; hold on; plot(OM(i+16,:),'color',cmap(i,:)/255); axis tight; title('PSP Inhibitory');
       subplot 414; hold on; plot(OM(i,:)+OM(i+8,:)+OM(i+16,:),'color',cmap(i,:)/255); axis tight; title('EEG = P+E+I');
    end 
    legend({'0','22.5','45','67.5','90','112.5','145','157.5'},...
        'Location','BestOutSide');
    
    % clear all
    clear; clearvars -global
end

function [OM] = VC_ode()
%% ODE solver using Euler method with constant stepsize 
global te ti he hi e0 r0 v0 nmass sz tend

% The connectivity and the weighting matrix
%disp('Loading the connectivity matrix...');
[W] = localConnectivity(); 

% Additional variables
ke = he/te; ki = hi/ti;
K  = [repmat(ke,2*nmass/3,1); repmat(ki,nmass/3,1)];
T1 = [repmat(-2/te,2*nmass/3,1); repmat(-2/ti,nmass/3,1)];
T2 = [repmat(-1/te^2,2*nmass/3,1); repmat(-1/ti^2,nmass/3,1)];

% Status counter 
trunk = 200; trunkidx = 1; 
n_trunk = ceil(tend/trunk/sz);

%%
%% Solver main loop
%%
OM = [];  
%Y0 = [zeros(nmass,1); zeros(nmass,1)];% Initial condition
%Y0(1) = -70; Y0(nmass+1) = -70;
Y0 = rand(2*nmass,1);
% Gradient
dv = zeros(nmass,1);
du = dv;

% Waitbar status
fprintf('Processing (total %d):    \t',n_trunk);

for idx = 1:length(0:sz:tend)      
    % Update half of Y0   
    u = Y0(1:nmass); 
    v = Y0(nmass+1:2*nmass);
   
    % Sigmoid function of membrane potentials    
    S = 2*e0 ./ (1+exp(r0*(v0-u))); 
    
    % Output
    S = W*S;  
       
    % ODEs 
    du = v;
    %dv(1:2*end/3)     = -2/te * v(1:2*end/3)     - 1/te^2 * u(1:2*end/3)     + S(1:2*end/3);
    %dv(2*end/3+1:end) = -2/ti * v(2*end/3+1:end) - 1/ti^2 * u(2*end/3+1:end) + S(2*end/3+1:end);          
    dv = T1.*v + T2.*u + K.*S;          
    
    % no noise
    % dY = [du; dv]
    Y1 = Y0 + sz * [du; dv];   
 
    % Update Y
    Y0 = Y1;
       
    % Save result in each step   
    OM = [OM u];
      
    if ~mod(idx,trunk)
       % display the informations
       % Waitbar status
       fprintf('\b\b\b\b%04d',trunkidx);
       trunkidx = trunkidx + 1;
    end

end

%% End of solver loops
fprintf('\n');    
end

function [W] = localConnectivity()
global p
    
%% The weighting matrix    
    c1s = p(7); c1n = p(7); c1l = p(9);
    c2s = p(10); c2n = p(11); c2l = p(12);
    
    wtMp = [c1s c1l c1l c1l c1l c1l c1l c1n;
            c1n c1s c1l c1l c1l c1l c1l c1l;
            c1l c1n c1s c1l c1l c1l c1l c1l;
            c1l c1l c1n c1s c1l c1l c1l c1l;
            c1l c1l c1l c1n c1s c1l c1l c1l;
            c1l c1l c1l c1l c1n c1s c1l c1l;
            c1l c1l c1l c1l c1l c1n c1s c1l;
            c1l c1l c1l c1l c1l c1l c1n c1s ];
        
     wtMi = [c2s c2l c2l c2l c2l c2l c2l c2n;
            c2n c2s c2l c2l c2l c2l c2l c2l;
            c2l c2n c2s c2l c2l c2l c2l c2l;
            c2l c2l c2n c2s c2l c2l c2l c2l;
            c2l c2l c2l c2n c2s c2l c2l c2l;
            c2l c2l c2l c2l c2n c2s c2l c2l;
            c2l c2l c2l c2l c2l c2n c2s c2l;
            c2l c2l c2l c2l c2l c2l c2n c2s ];  
        
%% The connectivity
    % The preferred angles
    angle_min = 0;
    angle_max = 180;
    d_angle   = 22.5;
    angle     = angle_min:d_angle:angle_max-d_angle;
    
    % Selectivity
    sel_mat = 1650.36; 
    
    % Orientation differences
    orientation_difference = [];
    for i = 1:8
        tmp = [];
       for j = 1:8
           tmp = [tmp abs(angle(i)-angle(j))];
       end
       orientation_difference = [orientation_difference; tmp];
    end
    orientation_difference(orientation_difference==157.5) = 22.5;
    orientation_difference(orientation_difference==135)   = 45;
    orientation_difference(orientation_difference==112.5) = 67.5;    
    
    % There is only local connection between the masses
    diffAngleLocal = orientation_difference;
    diffAngleLocal(diffAngleLocal>40) = NaN;
    %Wpp_local = sel_mat*0.01.*exp(-0.0072*diffAngleLocal); Wpp_local(isnan(Wpp_local)) = 0; 
    Wpp_local = sel_mat*0.002.*exp(-0.0072*diffAngleLocal); Wpp_local(isnan(Wpp_local)) = 0; 

    % P->eN = local + distal
    Wpe_local = Wpp_local;
    % P->iN = local + distal
    Wpi_local = Wpp_local;
    % Wp = [Wpp Wpe Wpi];
    Wp = [wtMp.*Wpp_local wtMp.*Wpe_local wtMi.*Wpi_local];

    % Excitatory interneuron 
    % eN->P = local
    Wep = Wpp_local;
    % eN->eN = local
    Wee = Wpp_local;
    % eN->iN = local
    Wei = Wpp_local;
    % We = [Wep Wee Wei];
    We = [wtMp.*Wep wtMp.*Wee wtMi.*Wei];

    % Inhibitory interneuron 
    % iN->P = local
    diffAngleLocal = orientation_difference;
    diffAngleLocal(diffAngleLocal>60) = NaN;

    %Wip_local = sel_mat*0.08.*exp(-0.038*diffAngleLocal); 
    Wip_local = sel_mat*5.*exp(-0.038*diffAngleLocal);
    Wip_local(isnan(Wip_local)) = 0; 
    Wip = Wip_local; 

    % iN-eN = local
    Wie_local = Wip_local;

    % iN-iN = local
    Wii_local = sel_mat*0.04.*exp(-0.038*diffAngleLocal); 
    Wii_local(isnan(Wii_local)) = 0; 

    Wi = [wtMp.*Wip wtMp.*Wie_local wtMi.*Wii_local];         
   
    W = [Wp; We; Wi];
    W = W';
end

function plotFigure(figname,Y)
    % Plot a time serie o
    Y1 = Y(:,end-ceil(4*end/5):end); 
    Y1 = Y1 - repmat(mean(Y1,2),1,size(Y1,2));

    % Color
    col = [1.0000  0  0; 0  0 1.0000;0  1.0000  0;0  1.0000  1.0000; 0.4784 0.0627 0.8941;
    0.6000  0.2000  0; 0    0.4980  0;  1.0000  0    1.0000];

    for i = 1:8 
          subplot(2,1,1); hold on; plot(Y1(i,1:end),'Color',col(i,:)); axis tight; title('Time series');
          subplot(2,1,2); hold on; plot(Y1(i,end-ceil(1*end/5):end),'Color',col(i,:));
          axis tight;  title('Last 1/5 of time series');   
          set(gca,'XTick',[1:length(ceil(size(Y1,1)/5))],'');
    end 
    
    % save figure;
    saveas(gcf,[figname '.png'],'png'); 
   % End
end

